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Strong numerical evidence is presented for the existence of a continuous family of time-periodic 
solutions with "weak" spatial localization of the spherically symmetric non-linear Klein-Gordon 
equation in 3-f 1 dimensions. These solutions are "weakly" localized in space in that they have slowly 
decaying oscillatory tails and can be interpreted as localized standing waves (quasi-breathers). By 
a detailed analysis of long-lived metastable states (oscillons) formed during the time evolution it is 
demonstrated that the oscillon states can be quantitatively described by the weakly localized quasi- 
breathers. It is found that the quasi-breathers and their oscillon counterparts exist for a whole 
\^ • continuum of frequencies. 
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I. INTRODUCTION 



Nonlinear wave equations (NLWE) lie at the heart of many fields in physics including hydrodynamics, classical 
integrable systems, Field Theories, etc. Let us give a prototype class of NLWE for a real scalar field, in 1 + n 
T— I . dimensional space-time: 
> ■ 

--| + A(")0 = F(0), where A(") ^ ^ , (1) 

a^ ■ 

' where the real function, F{(f)), defining the theory is given for example as = 0(0^ — 1) for the "canonical" 0^ 

model, and F{(f>) = sin((/)) for the Sine-Gordon (SG) model. 

A particularly interesting class of solutions of NLWE is the class of nonsingular ones exhibiting spatial localization. 
<^ , Such spatially localized solutions have finite energy and can correspond to static particle-like objects or to various 
"*7^ ■ traveling waves. In Field Theory localized static solutions have been quite intensively studied in a great number of 

models in various space-time dimensions, see e.g. the recent book of Sutcliffe and Manton [1]. 
^ Spatially localized solutions with a non-trivial time dependence (i.e. not simply in uniform motion) of NLWE are 

much harder to find. In fact the simple qualitative argument stating that "anything that can radiate does radiate", 
indicates that the time evolution of well localized Cauchy data leads generically to either a stable static solution plus 
radiation fields, or the originally localized fields disperse completely due to the continuous loss of energy through 
radiation. Since Derrick-type scaling arguments exclude the existence of localized static solutions in scalar field 
■ theories given by Eq. (1) in more than two spatial dimensions, for n > 2 one does not expect to find localized states 
at all at the end of a time evolution. Simplifying somewhat the above, one could say that in the absence of static 
localized solutions, localized initial data cannot stay forever localized. In fact for localized initial data there is a time 
scale which is the crossing time, Tc, (the time it takes for a wave propagating with characteristic speed to cross the 
localized region) and a priori one would expect the rapid dispersion of initially localized Cauchy data within a few 
units of Tc. One of the rare examples of a time-periodic solution, staying localized forever is the famous "breather" 
in the 1 + 1 dimensional Sine-Gordon (SG) model. 

In the 1 + 1 dimensional "canonical" (with a double well potential) (/)^ theory the pioneering work, based on 
perturbation theory by Dashen, Hasslacher and Neveu [2] indicated the possible existence of breather-like solutions. 
A completely independent numerical study by Kudryavtsev [3] has also indicated that suitable initial data evolve 
into breather- type states. These results stimulated a number of investigations about the possible existence of non- 
radiative solutions in the 1 + 1 dimensional 0* model. After a long history Segur and Kruskal [4] and Vuillermot 
[5] have finally established that in spite of the above mentioned perturbative and numerical indications time-periodic 
spatially localized finite energy solutions (breathers) do not exist. Even if genuine localized breathers in the 1 + 1 
dimensional 0^ model are absent, in view of the perturbative and numerical evidence for the existence of long living 
states "close" to genuine breathers, it is a natural question how to describe them. Boyd has made a detailed study 
of time periodic solutions which are only weakly localized in space (i.e. the field possesses a slowly decreasing 
oscillatory tail) but as long as the amplitude of the oscillatory "wings" are small they still have a well defined core. 
Boyd has dubbed such solutions "nanopterons" (small wings), we refer to his book for a detailed review [6]. 
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Quite interestingly Bogoluvskii and Makhankov have found numerical evidence for the existence of spatially localized 
breather-type states in the spherically symmetric sector of </)^ theory in 3+1 dimensions [7]. These breather- like objects 
observed during the time evolution of some initial data are called nowadays "oscillons" . Most of the observed oscillon 
states are unstable having only a finite lifetime. They lose their energy by radiating it (slowly) to infinity. More 
recent investigations started by Gleiser [8, 9] have revealed that oscillons do form in a fairly large class of scalar 
theories in various spatial dimensions via the collapse of field configurations (initial data) that interpolate between 
two vacuum states of a double well potential. Such a spherically symmetric configuration corresponds to a bubble, 
where the interpolating region is the bubble wall that separates the two vacuum states at some characteristic radius. 
These works have led to a renewed interest in the subject. 

It has been found that oscillons have extremely long lifetimes which is already quite remarkable and makes them 
of quite some interest. These long living oscillon states seem to occur in a rather generic way in various field theories 
involving scalar fields in even higher dimensional space-times [10] and according to [11] oscillons are also present in 
the nonabelian SU(2) bosonic sector of the standard model of elcctroweak interactions at least for certain values of 
the pertinent couplings. Such oscillons might have important effects on the inflationary scenario [12] as they could 
form in large numbers retaining a considerable amount of energy. In a recent study [13] of a 1 -|- 1 dimensional scalar 
theory on an expanding background exhibited very long oscillon lifetimes, while in Ref. [14], [15] extremely long living 
oscillons have been exhibited in a 1 + 2 dimensional sine-Gordon model. 

The sophisticated numerical simulations by Honda and Choptuik of the Cauchy problem for spherically symmetric 
configurations in the theory in 3 + 1 dimensions [16], [17], have revealed some interesting new features of the 
oscillons. In particular in Ref. [16] it has been found that by a suitable fine-tuning of the initial data the lifetime of 
the oscillons could be increased seemingly indefinitely, and it has been conjectured that actually an infinitely long 
lived, i.e. non-radiative, spatially localized solution exists. Furthermore the existence of such a solution would even 
provide the explanation of the "raison d'etre" and of the observed genericness of long lived oscillons. The eventual 
existence of a non-radiative breather in this simple and "generic" (p* model in 3 + 1 dimensions would be clearly 
of quite some importance. It should be noted, that in quite a few spatially discrete models, localized time-periodic 
"discrete breathers" have been shown to exist and they are being intensively studied [18]. 

One of the motivation of our work has been to clarify if non-radiative breathers indeed exist and to find them directly 
by studying time-periodic solutions of the NLWE. Our numerical results led us to conclude that no localized (with 
finite energy) time-periodic solutions exist in the 0^ model. On the other hand this study has led us to understand the 
oscillon phenomenon better and we present a simple but quantitatively correct scenario explaining some important 
properties of the oscillons (such as their existence and their long lifetimes). Our scenario is based on the existence 
of a special class of time periodic solutions which are weakly localized in space. Such solutions (which have infinite 
energy) will be referred to as quasi-breathers (QB). 

In the present paper we make a detailed study of oscillons in the (already much studied) cj)'^ model in 3 -I- 1 
dimensions. Using a previously developed and well tested time evolution code where space is compactified, thereby 
avoiding the problem of artificial boundaries [19], we compute some long-time evolutions of Gaussian- type initial 
data. We observe that long living (6000 — 7000 in natural units) oscillon states are formed from generic initial data. 
These oscillons radiate slowly their energy and for short (as compared to their total lifetime) timescales they can be 
characterized by a typical frequency, uj. This frequency increases slowly during the lifetime of the oscillon and when 
u) reaches a critical value ~ 1-365 there is a rapid decay. By fine-tuning the initial data one can achieve that the 
oscillon state instead of rapidly decaying evolves into a near time-periodic state, whose frequency is nearly constant 
in time and 1.365 < lu < 1.412. The existence of such near-periodic states (referred to as resonant oscillons) has been 
already reported in Ref. [16], and we also observe an increase of the lifetime of these states without any apparent 
limit by fine-tuning the parameters of the initial data to more and more significant digits. There are, however, also 
some discrepancies. We find that clearly distinct near-periodic states for various values of the pulsation frequency 
1.365 < UJ < 1.412 exist. According to our results there is little doubt that for any value of uj in this range a 
corresponding near-periodic oscillon state exists. Our data clearly shows that the near-periodic states also radiate, 
although very weakly. The radiation becomes weaker and weaker as w ^ v^. 

On the other hand we have implemented a multi domain spectral method in order to find directly stationary, time 
periodic solutions of the NLWE and compare them to the long living oscillon states obtained from the time evolution. 
This makes it possible to attack the problem of finding directly the putative time periodic breather(s). (One could 
easily generalize our method for the quasi-periodic case.) We find that there is a large family of time periodic solutions 
which are only weakly localized in space, in that they have a well defined core, and an oscillatory tail decreasing as 
DC 1/r. We single out a special family among them by minimizing the amplitude of their oscillatory tail which definition 
comes close to minimizing the energy density of the oscillatory tail. These solutions are the closest to a breather and 
for that reason we call them quasi-breathers. They seem to exist for any frequency, < uj < \/2, although in this 
paper we exhibit QBs only with 1.30 < lo < \f2 . The amplitude of the oscillatory tail of the QBs becomes arbitrarily 
small as the frequency approaches the continuum threshold defined by the mass of the field, w — > \/2. Our numerical 
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evidence speaks clearly against the possible existence of a truly localized, breather-like solution periodic in time, for 
the frequency range 1.30 < lu < 1.412 contrary to the claims of Ref. [16]. In view of these conflicting numerical findings 
it is now highly desirable to try to find an analytical proof or disproof of the existence of a localized non-radiative 
solution to settle this issue. We do not expect the situation being qualitatively different from the 1-1-1 dimensional 
case, and although the proofs of Refs. [4, 5] are not applicable for the 3-1-1 dimensional case, we see no reason that 
their negative conclusion would be altered. 

More importantly, we believe to have made a step towards understanding the mechanism behind the existence of such 
long-living oscillonic states without the need to invoke genuine breather solutions, which even if they would exist would 
be clearly non-generic, while the QBs seem to be generic. The total energy of the quasi-brcathers is divergent, due to 
the lack of sufficient (exponential) spatial localization, hence they are not of direct physical relevance. Nevertheless 
a careful numerical study shows that the oscillons produced during the time evolution of some suitable Cauchy data 
arc quantitatively very well described by the quasi-breathers. By comparing the Fourier decomposition of an oscillon 
state at some instant, t characterized by a frequency, Lo{t), obtained during the time evolution with that of the 
corresponding QB, we have obtained convincing evidence that the localized part of the oscillon corresponds to the 
core of the QB of frequency u!{t). What is more, the oscillatory tail of the QB describes very well the standing wave 
part of the oscillon. Our oscillon scenario is based on this analysis and leads us to propose that any oscillon contains 
the core and a significant part of the oscillatory tail of the corresponding QB. The time evolution of an oscillon can 
be approximatively described as an adiabatic evolution through a sequence of QBs with a slowly changing frequency 
u>{t). This oscillon scenario is based on the existence and of the genericness of the QBs. 

This paper is organized as follows. In Sec. II we study the time evolution of localized, Gaussian-type initial data in 
(/>^ theory and investigate some important aspects of the oscillon solutions. In Sec. Ill we present an infinite system 
of coupled ordinary differential equations (ODE's) obtained by the Fourier-mode decomposition of the NLWE Eq. 
(3) and discuss some of its properties. Section IV is devoted to the description of the spectral methods used to solve 
this system. In particular, we carefully deal with the asymptotic behavior of the Fourier-modes. Various convergence 
tests are exhibited. The quasi-breathers are discussed in Sec. V, the results on our oscillon scenario are discussed in 
Sec. VI and conclusions are drawn in Sec. VII. 



where 4> is a real scalar field, dt = d/dt, di = djdx^ and i = 1,2,3. In this paper we shall restrict ourselves to 
spherically symmetric field configurations, when the corresponding NLWE is given by 



II. TIME EVOLUTION 



A. The nonlinecir wave equation of the 0* theory 



We consider the following theory in 1 -I- 3 dimensions whose action can be written as 




(2) 




where A 




(3) 



The energy corresponding to the action (2) can be written as 




(4) 







where ^ denotes the energy density 



(5) 



It is easy to see that the finiteness of the total energy is guaranteed by </> ^ ±1 -|- ^^(r ^/^) as r ^ oo. 



4 



B. Numerical techniques 



We briefly outline here the main ideas for the implementation of our evolution code to solve numerically the Cauchy 
problem for Eq. (3). Assuming ^ — > 1 as r — > oo we introduce the new field, (j) as 

^(,.) = ^-l. (6) 
r 

Then the NLWE Eq. (3) takes the form: 

"2 (dU - dU) = ${4>-r)[4>- 2r) . (7) 



The next step is to compactify in the spacelike directions by a suitable coordinate transformation of r. This way 
we can guarantee that our computational grid, associated with a finite-difference scheme, covers the entire physical 
spacetime, at least in principle. Specifically, a new radial coordinate, R, is introduced in the following way: 

where K is an arbitrary positive constant. In the new radial coordinate, R, the entire Minkowski spacetime is covered 
by the coordinate domain < i? < 1 while spacelike infinity is represented by the 'hypersurface' R= 1. The R = const 
'lines' represent world-lines of 'static observers', i.e. integral curves of the vector field {d/dt)"". In the compactified 
representation the field equation, (7), reads as 



i?2 



where Cl is given by 



n = -(i-R'). (10) 



We remark that the spacelike compactification used here is a simplified variant of the conformal transformation 
used in [19]. There instead of the t = const, hypersurfaces the initial data are specified on hyperboloids. Furthermore 
Minkowski spacetime is compactified mapping null infinity to finite coordinate values. Since in the present case the 
scalar field, (j), is massive, i.e. never reaches null infinity the hyperboloidal compactification is not essential. 

In order to obtain a system of first order equations we introduce the independent variables = dfcf) and t^/j = Or^. 
Then equation (9) can be rewritten as 



K2(l + i?2)2-«-^« K(l + i?2)3 jyj^ 

which together with the integrability condition 9*0/? = dji^'t and the defining equation dt^ = (f)t form a strongly 
hyperbolic system of first order differential equations for the three variables 0, cf)t and (f)R (sec e.g. [20]). The initial 
value problem for such a first order system is known to be well-posed [21]. Note that the relation dB,<t> = is 
preserved by the evolution equations, and therefore it corresponds to a constraint equation. 

In order to solve the initial value problem for equation (11) we discretize the independent variables t and R. A 
simple uniform grid with steps Ai and Ai? is introduced. Spatial derivatives are calculated by symmetric fourth order 
stencils. Time integration is done using the 'method of lines' in a fourth order Runge-Kutta scheme, following the 
recipes proposed by Gustafsson et al [21]. A dissipative term proportional to the sixth derivative of the field is added 
in order to stabilize the evolution. Since this dissipative term is also chosen to be proportional to (Ai?)^, it does not 
reduce the order of the applied numerical method, in other words, its influence is decreased by the increase of the 
used resolution. The numerical methods and the actual numerical code we use for calculating time evolution in this 
paper are also based on those developed in [19]. 

A few non-physical grid points are introduced for both negative radii i? < and for the region "beyond infinity" 
R> 1. Instead of calculating the time evolution of the i? < points, the symmetry property of (j) about the origin 
i? = is used to set the function values at each time step. Similarly, (!>, being a massive field, decays exponentially 
towards infinity, consequently all the field values (j)^ (j}i and 4>j^ are set to zero for i? > 1 during the entire evolution. 
This takes care of the spacelike infinity = in equation (11). Therefore it is possible to use symmetric stencils 
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exclusively when calculating spatial derivatives. The grid point at the origin R = needs special treatment since 
the last term as it stands on the right hand side of equation (11) is apparently singular. However, this term, when 
evaluated in terms of the original (non-singular) field value (j) have zero limit value at i? = 0. 

Although compactifying in spatial direction restricts the coordinate i? to a finite domain, grid points in our numerical 
representation get separated by larger and larger physical distances as we approach R = 1. This far region, where 
shells of outgoing radiation cannot be represented properly, moves out to higher and higher physical radii as AR 
decreases. Numerical simulations with increasing number of grid points demonstrate that wave packets of outgoing 
massive fields get absorbed in the transitional region without getting reflected back into the inner domain. In this 
way our numerical simulations still give a good description of the field behavior precisely in the central region for very 
long time periods. The simple but physically non-uniform grid together with the help of the dissipation term appears 
to absorb outgoing radiation in a similar way to the explicit adiabatic dumping term method applied by Gleiser and 
Sornborger [22]. Furthermore, because of the very low inward light velocity in the asymptotic region i? w 1 the field 
behavior in the central region is correctly given for a long time period even after the appearance of numerical errors 
at i? w 1. 

Simulating time evolution of oscillons up to their typical maximal lifetime of f = 7000 using spatial resolution of 
2^^ points took usually a week on personal computers. However, because of the need of several runs when fine-tuning 
parameters in the initial data, we mostly used typical resolutions of 2^^ spatial points. The parameter k in the 
coordinate transformation (8) was set to k = 0.05 in our simulations in order to concentrate approximately the same 
number of grid points to the central oscillon region and to the far away region where the massive fields form high 
frequency expanding shells. A Courant factor of ^ = 1 turned out to be appropriate to obtain stable simulations 
with our choice of k. 

The convergence tests confirmed that our code does provide a fourth order representation of the selected evolution 
equations. Moreover, we monitored the energy conservation and the preservation of the constraint equation dfj4> = 4>ri.- 
Most importantly, we compared the field values which can be deduced by making use of the Green's function and by 
the adaptation of our particular numerical code to the case of massive linear Klein-Gordon fields [23]. The coincidence 
between the values in the central region provided by these two independent methods for long time evolutions (i « 10^ 
measured in mass units) made it apparent that the phenomena described below should be considered as true physical 
properties of the investigated non-linear field configurations. 

C. Oscillons 

Following Refs. [9, 16] we start with the following Gaussian-type initial data: 

</>|t=o = 0OO + (0c - (f>oo) ■ exp(-rVro) , dt(p\t^„ = , (12) 

with (j)c and being the field values at the center r — and at infinity r = oo while tq is the characteristic size of 
the bubble at which the field values interpolates between (pc and (poo- By fixing (poo = —1 as in [9, 16] but varying 
ro and 4>c, Eq. (12) provides a two-parameter family of smooth and suitably localized initial data. For a large open 
subset of the possible initial parameters (pc and ro, after a short transitional period the field evolves into a long living 
localized nearly periodic state, named oscillon by Copeland, Gleiser and Miiller. Although these configurations five 
much longer than the dynamical time scale expected from the linearized version of the problem (i.e. light crossing 
time) , their lifetime is clearly not infinite. The energy of these oscillating states is slowly but definitely decreasing in 
time, and after a certain time period they quickly disintegrate. For the time dependence of the energy in a compact 
region see Fig. 3 of [9]. We illustrate on Fig. 1 two such oscillon states with rather different lifetimes. 

The parameter dependence of the lifetime is illustrated on Figs. 6, 7 and 8 of [9]. We note that since the final 
decaying period is relatively short, furthermore its time dependence is almost the same for each particular choice of 
initial data, the lifetime plots are quite insensitive of the precise definition of how one measures the lifetime of a given 
configuration. In our calculations the lifetime r was defined by observing when the value of the oscillating field at the 
center r = falls (and remains) below a certain prescribed value (e.g. (p = —0.95). 

Already Copeland, Gleiser and Miiller have noticed that a delicate fine structure appears in the lifetime plot. This 
peculiar dependence on the precise value of the initial parameter ro has motivated the detailed investigation of Honda 
and Choptuik (see Figs. 4 and 5 of Ref. [16]). The calculations have shown that fixing the value (pc the lifetime 
increases without any apparent upper limit when the parameter ro approaches some element of a large set of discrete 
resonance values r*. For example, in the case (pc = 1 Honda and Choptuik have found 125 peaks on the lifetime plot 
between 2 < ro < 5. These fine-tuned oscillons at a later stage during their evolution develop into a state which is 
very close to a periodic (non-radiating) one. On Fig. 2 we show the field value at the center for typical oscillons 
close to a chosen peak. The seemingly stable almost periodic stages of the fine-tuned oscillon configurations will be 
referred as near-periodic states. The closer the initial parameters are to the critical value the longer the lifetime of the 
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FIG. 1: The upper and lower envelope of the oscillations of the field 4>{r = 0) is shown for the evolution of Gaussian initial 
data of the type (12) with two different sets of initial parameters. It is important to note that the field oscillates with a 
non-constant period 4.6 < T < 5. For generic oscillons the period T shows a decreasing tendency towards T « 4.6. 
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FIG. 2: Upper envelope of the field value for different oscillon states close to the 27th peak of the lifetime curve belonging to 
4>c = 1. These states are characterized by initial value ro which is near rjy « 2.39297. The near-periodic state is approximately 
at 2100 <t< 2800 for the fine-tuned oscillons in the lower plot. 
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near-periodic state becomes. Although the period T of the underlying high frequency oscillations remains constant 
to a very good approximation during a chosen near-periodic state, different near-periodic states, corresponding to 
various peaks on the lifetime curve, oscillate with clearly distinct periods in the range 4.446 < T < 4.556. Generic 
oscillon states and even the not near-periodic initial part of fine-tuned oscillons pulsate with longer period, in the 
range 4.6 < T < 5. 

Generic oscillons with initial data between two neighboring resonance values r* < ro < r*^^ may have already quite 
a long lifetime without developing into a near-periodic state. Plotting the field value at the center r = one can 
see a low frequency modulation of the amplitude (called shape mode by Honda and Choptuik) of the high frequency 
basic oscillating mode (see upper plot of Fig. 2). Oscillons between the next two resonance values r*j^i < < r*_|_2 
are distinguished from those in the previous interval by the fact that they possess exactly one more or one less peak 
associated with these low frequency oscillations on the envelope of the field value <^(t, 0) before they disperse. Longer 
living supercritical states arise when one closely approaches a critical value from one of the possible two directions. 
Then the last low frequency modulation peak gets shifted out to a later and later time as one goes towards the 
resonance, making room for a near-periodic state between the last two modulations. Close to the resonance, but on 
the other side of it, the same long living near-periodic state appears, now called subcritical, with the only difference 
that at the end the field disperses without forming a last low frequency modulation peak. 

D. Results 

It was shown in Ref. [16] that close to a resonance the oscillon lifetime r obeys a scaling law 

Tr^'yln\ro-r*o\+6 , (13) 

where the scaling exponent 7 has specific values for each resonance, while the constant 5 is also different for subcritical 
and supercritical states. Although the lifetime appears to increase without any limit by fine-tuning the initial pa- 
rameter, in practice it is very difficult to achieve very long lifetimes because one cannot represent numbers very close 
to the resonance value due to the limitation implied by the applied machine precision. Achieving longer lifetime for 
the near-periodic state is possible by using high precision arithmetics, although then there is a considerable increase 
of the computation time which limits the applicability of this approach. Using "long double" variables on an SGI 
machine, or applying the software "doubledoublc" [24] on a personal computer, we could calculate with twice as many 
significant digits than standard double precision computer variables can represent (i.e. 32 instead of 16). This way we 
could improve the fine-tuning and double the observed lifetime of near-periodic states. On Fig. 3 we plot the scaling 
law of the lifetime for oscillons near three different resonances. Instead of choosing three resonances with 0c = 1, we 
calculated states close to the first peak on the lifetime curves corresponding to three different values of (pc- The reason 
for this was that the normal, slowly but evidently decaying, oscillon state is the shortest near the first peak (i.e. no 
modulation on the contour curve), thereby we could concentrate computational resources on the near-periodic state. 

In order to clarify what we mean by fine-tuning the parameter ro to 16 (or 32) digits and what is the actual error 
of the quantities, we present a table on the precise location of the first peak for cpc = ^ when performing fine-timing 
with five different numerical resolutions. For each spatial resolution we could achieve approximately the same lifetime 
of approximately r = 1100 when the parameter ro approximated a resolution dependent value to 16 digits. The 
convergence of the data indicates that the actual position of the peak is at Tq = 2.281370594 with an error of 10~^. 
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TABLE I: The position Tq of the first peak for 0e = 1 using various resolutions. The number of spatial grid points used 
for a specific fine-tuning is Ui = 2\ Tlie error is estimated as 5i = Itq*-'^ — Tq^^^Y The convergence factor is defined as 

Our numerical simulations clearly show that the different resonance peaks on the lifetime plot correspond to different 
near-periodic states. In fact, a one-parameter family of distinct near-periodic states appears to exist. This statement is 
in marked contrast with the claim of Honda and Choptuik in [16], where in Sec. III-A they claim that the oscillation 
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FIG. 3: The oscillon lifetime, r, versus — In |ro — ro\ is shown for states close to three different resonances. The lower lines 
plotted with a given line type represent subcritical, ro < rg, states while the upper lines with the same line type correspond to 
supercritical , ro > ro, solutions. 



period is almost the same for all oscillons and is roughly T ~ 4.6, which corresponds to a pulsation frequency of 
CO ~ 1.366. Prom our numerical analysis we can clearly see that the period of the oscillation depends on the resonance 

considered. On Fig. 4 the upper envelope of the oscillating field value at the origin is shown as function of time for 
oscillon states close to the same three resonances as on Fig. 3. On Fig. 5 the time dependence of the frequency at the 
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FIG. 4: Upper envelopes of the oscillations of the field value at r = for three pairs of sub and supercritical oscillon states 
close to the first resonances corresponding to three given 4>c initial parameters. The fine-tuning was performed using 32 digit 
arithmetics. 



origin is shown for the same three states during the time period where the oscillation is almost periodic (i.e. during 
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the near-periodic states). The frequency w of the oscillations along with its time dependence at some radius r = r has 
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FIG. 5: The pulsation frequency, u) in function of the time for the three pairs of oscillons shown on the previous figure during 
their corresponding near-periodic states. 

been determined from our numerical results by minimizing the following integral for the oscillation period T = 27r/a; 



rt+to 

/ [<j>{t,f)-<P{t + T,f)X 

Jt-to 



dt 



(14) 



using some suitably chosen integration interval determined by ^o- This procedure, applying polynomial interpolation, 
yields significantly more precise frequency values than the direct use of the Fast Fourier Transform method, especially 
when the time step by which our evolution code outputs data is not extremely short. Another advantage of our 
procedure is it being much less sensitive to the particular choice of the sampling interval (in this case 2to) than the 
FFT algorithm. The relatively small value of to (for example to = 10) makes it possible to monitor relatively sharp 
changes in the time dependence of the frequency. 

All three near-periodic states show a low frequency change of u> with a decaying amplitude, with the lowest frequency 
modulation belonging to the state with the highest ui. The amplitudes on Fig. 4 would show a similarly decaying slight 
modulation if we would plot them individually in time intervals where the oscillations are almost time-periodic. These 
observations suggest that near-periodic states of frequency to also contain a superposition of states with frequencies 
ui ± A^j where A^/lo ^ 1. In the next sections we shall see that at least the core part of an oscillon of frequency u) 
can be extremely well described by a weakly localized quasi-breather of the same frequency. 

We emphasize, that the three chosen near-periodic states are typical, i.e. the near-periodic states of all fine-tuned 
oscillons (including all peaks belonging to (j)c = 1) are qualitatively and quantitatively very similar to them. Actually, 
the frequency of all calculated near-periodic states fell into the interval [1.379, 1.413] spanned by the frequency of the 
two chosen first peaks belonging to 0c = —0.4 and (f)c = —0.7. 
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On the third plot of Fig. 5 we can also sec a slow but steady decrease of the frequency oj. For the other two states, 
with frequencies closer to \/2 no such behavior is apparent in the time interval we could simulate. This slow decrease 
is also manifested in the energy of the configuration. On Fig. 6 we plot the time dependence of the energy contained 
inside spheres with three subsequent radii. The decrease of the energy indicates that the near-periodic state slowly 
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FIG. 6: The time dependence of the energy contained in spheres of radii r — 19.5, r — 40.1 and r = 107.1 for the (j)c = —0.4 
first peak. The thickness of the curve indicates that the high frequency oscillation is still not negligible at these high radii. The 
different endings of the curves correspond to sub and supercritical states. 

loses its energy, consequently it cannot be exactly time-periodic. Looking at the energy in the sphere at r = 40.1, we 
can see that the decrease of energy from t = 1000 to t = 2000 is AE = 0.0041. Taking into account that the total 
energy is -E = 21.60, we can give a naive estimate on the lifetime by calculating when the energy would decrease into 
its half value at this rate, getting Tg = 2.6 • 10^. We expect that all near-periodic states radiate, but this radiation is 
becoming considerably weaker as the frequency get closer to \/2. This expectation is confirmed by the direct analysis 
of periodic solutions in the next sections. 

On Fig. 6 the curve corresponding to the largest radius (r = 107.1) indicates how slowly the energy moves outwards 
becaiise of the massive character of the field. The different endings of the curves belonging to sub and supercritical 
states illustrate the two kinds of decay mechanism of near-periodic states. In the first, subcritical mechanism, the field 
quickly disperses, while energy moves essentially outwards only. In the second supercritical way, the energy of the 
field first collapses to a smaller region near the origin and then disperses to infinity. This resembles to the behavior 
of some unstable spherical shell, although no shell structure is visible on the density plots, being the highest energy 
density always at the center. The instability of the near-periodic states with two distinct kinds of decay mechanisms 
gives a qualitative explanation on why it is possible to reach long lifetimes by fine-tuning the initial parameters. 

E. Fourier decomposition of the evolution results 

Since during the time evolution the field (j) becomes approximately time-periodic for any longer living oscillon state, 
it is natural to look at the Fourier decomposition of the results provided by our evolution code. Since the Fast 
Fourier Transform algorithm is very sensitive to the size of the time step with which our evolution code writes out 
data, we use an alternative direct method which turns out to be significantly more precise in determining the basis 
frequency and the amplitude of the higher modes. As a first step we determine the oscillation period by locating two 
subsequent maxima, at instants ti and t2, of the field at the origin r = 0. Since ti and t2 fall in general between 
two conseqiient time slices written out by our evolution code, we approximate their position by fitting second order 
polynomials on the data. It is apparent from our results that for near-periodic states these maxima correspond to 
two consecutive moments of time symmetry not only at the center but also in a large region around the center to a 
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very good approximation. After determining the oscillation frequency lo = 27r/(t2 — ^i), we obtain the n-th Fourier 
mode at a radius r by calculating the following integral using the output 4>{t, r) of the evolution code 

4>n{r)= I {(j){t,r) + l)exp{inuit)dt . (15) 

We note that care must be taken for the first and last incomplete time steps when evaluating the integral (15). In 
the case of exact periodicity and time symmetry the imaginary part of this integral would be zero for all n. In order 
to check the deviation from time symmetry for various radii at the moments ti and ^2 we also evaluate the imaginary 
part of the integral for each n and verify whether it is small compared to the real part. 

As an example, on Fig. 7 we give the Fourier decomposition of the t^c = —0.7, ro = 9.38283 near-periodic state 
presented on Fig. 5. The integral (15) was calculated between two subsequent maximums of the field at the center 
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FIG. 7: Radial dependence of Fourier modes obtained by decomposition of a near-periodic state with frequency u) = 1.41203. 

at tl = 5420.22 and t2 = 5424.67, yielding an oscillation frequency of u = 1.41203. For the decomposition of a 
near-periodic state with a lower frequency see Fig. 9 of [16], where the frequency is a; = 1.38 (page 107. of [17]). 
From the two plots one can already see the general tendencies. As the basis frequency increases towards the upper 
limit the oscillon becomes wider, although with a decreasing amplitude. The infiuence of the higher modes (the 
relative amplitude ratio) is also getting smaller when the frequency grows. 



III. FOURIER DECOMPOSITION OF THE QUASI-BREATHERS 

This section is concerned with a direct search for time-periodic solutions of the NLWE Eq. (3) by means of a 
Fourier-mode decomposition of the scalar field cf). A simple parameter counting already indicates that an infinite 
number of parameters would be necessary in order to obtain an exponentially localized breather, so even if some do 
exist, it seems to be difficult to produce them directly. Therefore we look for weakly localized solutions of Eq. (3), 
and by analyzing their behavior in function of the free parameters we expect to obtain a clear indication about the 
existence of truly localized breathers. 
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A. Equations and their Asymptotic behaviors 

Since we seek periodic solutions of Eqs. (3) we make a Fourier mode decomposition of the form 

oo 

(i, r) = -1 + ^ = -1 + 00 (r) + ^ 0„ (r) cos {ncot) . (16) 

We have assumed here that by a suitable choice of the time origin the solutions are time-symmetric. Inserting this 
form in Eq. (3) gives rise to a system of coupled elliptic equations 

oo ^ oo 

[A - 2] (^0 = 00 (^0 - 3) + - (00 - 1) ^ 0m + ^ X] 't>m(l>p(l)q 5m,±p±q (17) 

m=l m,p,q—l 

- ^n] = 300 (00 - 2) 0„ + - (00 - 1) 0m0p 5„,±7n±p 

m,p=\ 

^ oc 

+ ^ ^ 0m0p09 f^n.itmibpibg 

where = (2— n^w^). In fact we have already regrouped together all the linear (first order) terms of the corresponding 
Fourier mode, 0„ to the left hand site of Eqs. (17). Wc will make a systematic search for localized solutions of the 
system (17). We shall not assume any a priori connection with the oscillons presented in the previous section for the 
possible frequencies, w. As we shall see below the eventual existence of a localized solution can only be possible if 
system (17) exhibits some very special properties. 

In order to look for regular solutions of Eqs. (17) we analyze the various asymptotic behaviors of the homogeneous 
solutions of the corresponding operator on the left hand side of Eqs. (17). 

• If > then there is one regular homogeneous solution, which falls off exponentially : exp [— |A„|r] /r. This 
function clearly decays sufBciently fast at infinity so that the energy would always be convergent. 

• A^ = is a degenerate case when the operator is the usual Laplacian. The non-singular (decaying) homogeneous 
solution is simply 1/r which, however, docs not decay sufficiently fast for the energy to be bounded. 

• If A^ < 0, both homogeneous solutions tend to zero at infinity, cos[|A„|r]/r and sin [|A„|r] /r. As for the 
previous case, these functions do not decrease fast enough to make the energy convergent. Apart from that, 
both functions arc regular at infinity and this implies that there is no uniqueness of the solution which is only 

defined up to a given phase (for more details, see Sec. IV C). 

Observe that no matter what the value of uj is, there always exists an n^^, such that aiJ modes, (pn n > tIoj will be 
of the third type (i.e. with A^ < 0). It is now easy to understand why in the generic case, one cannot expect to 

find exponentially localized ( with finite energy) solutions of Eqs. (17). To ensure localization one has to suppress 
all slowly decaying oscillatory modes which would normally require an infinite set of freely tunable parameters. The 
problem of finding a localized solution can be seen as a matching problem between the set of modes regular in a 
neighborhood of the origin, {0°}J^o ^^"^ modes with fast (exponential) decay for r — > oo, {05f }5^o- 

this case each mode 0*^ regular at the origin, has a single freely tunable parameter whereas none of the modes with 
fast fall off, 0^ for n > have any free parameters and therefore a sort of a miracle is needed that the two sets 
could be matched. This counting implies that while one can expect to find time-periodic solutions Eqs. (17) (even a 
whole family), these solutions have generically oscillatory tails for large values of r. This clearly reflects the argument 
"anything that can radiate docs radiate" transposed to the stationary case. On the other hand, it is not excluded 
that a localized solution might exist for very particular values of w, for which the oscillatory tails are absent. 

IV. NUMERICS FOR THE MODE DECOMPOSITION 
A. PDE solver 

We solve the system (17) using the LORENE library [25]. The basic features of LORENE are the use of spectral 

methods and multi-domain decomposition. The problem we are facing here being purely spherical, the fields are 
expanded on Chebyshev polynomials. The physical space is decomposed in various spherical domains. 
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Using such techniques, solving differential equations can be reduced to, in each domain, inverting a matrix on the 
coefficients space. Then, a linear combination of the particular solutions with the homogeneous ones is done, in order 
to impose regularity at the origin, appropriate boundary conditions and continuity of the overall solution. We refer 
the reader to [26] for more details on the algorithm, in the case of a Poisson equation. 

For this work, we have extended the operators presented in [26] and included the Helmholtz operators appearing 
in Eqs. (17), i.e. A — with both signs of A^. This is rather straightforward because, as for the Laplacian, they 
produce two homogeneous solutions. Only the case of infinity has to be treated differently as we will see in Sec. IV C. 

B. Description of the sources 

As wc have seen in Sec. Ill A, for all value of uj, there exists a value of n after which all the modes arc dominated by 
homogeneous solutions of the type sin (|A„| r + (fin) /r- Such functions are not easily dealt with by our solver. Indeed, 
in order to treat spatial infinity, LORENE usually uses a simple compactification by means of the variable u= 1/r in 
the exterior domain. In the past, this has enabled us to impose exact boundary conditions at infinity. When waves 
are present, it is well known [27, 28] that such a compactification is not practical. Indeed, no matter what scheme is 
used, there is always a point after which the distance between computational points (grid or collocation) is greater 
than the characteristic length of change of the wave, thus causing the scheme to fail. This issue is dealt with by 
imposing boundary conditions at a finite radius i?iim. for the slowly decaying oscillatory modes. The reader should 
refer to [29] where such methods have been applied to gravitational waves. For the exponentially decaying modes 
in Eqs. (17) we decide to keep in the right hand side (in the sotirces) only those terms that are dominated by the 
exponentially decaying homogeneous solutions, i.e., only those modes for which A^ > 0. As it will be seen later, the 
range of interesting pulsations is a; < \/2. Therefore the only two exponentially decaying modes arc; c6o and <J)i. In 
the Eqs. for 0o and (/>i for r > Rum., we set all the higher modes (n > 2) to zero (including terms of the type ^3(/>f ). 
So, we effectively solve for large values of the radius the following equations: 

[A-2]<^o = ^(</'o-l)<^i+0o(<^o-3) (18) 
[A - A^] = 3</)o ((/)o - 2) + ^<j)l . (19) 

This method yields solutions for (/)o and 0i which are correct for "intermediately large" values of r > -Riim. where the 
oscillatory and slowly decaying terms induced by the nonlincaritics do not dominate. It is clear that for sufficiently 
large values of r (j)Q and do not decay exponentially, since their behavior will be dominated by the slowly decaying 
oscillatory nonlinear source terms. We have carefully checked that changing the value of i?iim. does not influence the 
oscillatory modes, therefore we can conclude that the back reaction of the (/)o and is negligible on n > 2 for 

r > iJlim.- 

C. The operators 

When A^ > 0, the Helmholtz operator A — A^ admits two homogeneous solutions of which only one tends to zero 
at infinity. This situation is exactly the same as the one for the standard Laplace operator and all the techniques 
presented in [26] can be used. Once again, as we will be working for w < \/2, this happens only for and 4>\ and 
the associated sources have been given in the previous section IV B. 

The situation is quite different when dealing with the Helmholtz operator with A^ < 0. Apart from the compactifi- 
cation problem previously mentioned, wc have to note that now there are two homogeneous solutions that are regular 
at infinity. There is no reason to prefer one to the other. This means that there is no unique solution to our problem. 
Indeed, one can get a solution by doing the matching with any homogeneous solution of the type sin (|A„| r + (p„) /r, 
where can take any value in [0, 2ti[. So to summarize, for all modes such that A^ < we have to match, at a finite 
radius i^um. , the solution with a homogeneous one of the type : 

0n>2 (r > -Rlim) = AnS\n{\\n\r -\- ifin) jr. (20) 

Clearly we need additional conditions to fix the values of the phases <^„. In order to do so, let us recall that we 
are mainly interested in finite energy solutions. Such solutions should not contain any oscillatory behavior in 1/r 
at infinity, i.e. ell the coefficients A„ of such homogeneous solutions should be zero. One can hope to achieve that 
by searching, in the parameter space of the phases ((^2, <^n), the values that minimizes the absolute value of 



14 



cocfBcient of the first oscillatory homogeneous solution that appears : \A2\. This solution being "quite close" to a 
localized breather, for that reason we refer to it as a quasi-breather. 

Given the non-linearity of the problem, the location of the minimum can not be found analytically and one has to 
rely on a numerical search. We make use of the multi-dimensional minimizer provided by the GSL numerical library 
[30]. The algorithm is based on the simplex algorithm of Nelder and Mead [31]. We start by setting all the phases to 
7r/2 and iterate the procedure until the value of the minimum converges with a given threshold. Let us mention that 
the phases are searched in [0, 7r[, the rest of the interval, being described by changing the sign of the amplitudes. The 
simplex solver converges rapidly, the function \A2\ being very smooth, with no local extrema. 

So the final situation is the following : 

• For (f)o and (f)i the space is compactified and the sources in the exterior region are given by Eqs. (18) and (19). 

• For (j)n, n > 2, we only solve for r < and match the solution with a homogeneous solution 
An sin (|A„| r -|- 1^„) /r, the being determined by the simplex solver by minimizing |A2|. 

D. Avoiding the trivial solution 

The system (17) is solved by iteration but a problem arises from the fact that the trivial solution = is a 
solution of the equations (then 0o is a constant whose value is either 0, 1 or 2). No matter what the initial guess for 
the different modes is, the code always converges to a trivial static solution of this type. 

So, one needs to find a way to prevent this to happen. Given that the system (17) is coupled, it is sufficient to 
impose (f)\ (r = 0) ^ 0, in order to avoid the trivial solution. To do so, after each step of the iteration, we rescale (t)\ 
everywhere, by a factor a in order to impose that (^i (r = 0) has a certain value. In the general case, after convergence 
of the code, this scaling parameter is different from 1 meaning that the value of 4>i we found is no longer solution of 
the system (17). However for some values of u (being in an interval, see Sec. V), it is possible to find, exactly one 
value of (r = 0) such than the scaling factor is one. Thus, it is possible, at least for some values of u) to find the 
appropriate value of (f)i (r = 0) such that the obtained modes are non-trivial solutions of Eqs. (17). 

In practice, after convergence to a threshold level of typically 10~^ we switch on the convergence toward the true 
value of 01 (r = 0). With a technique already use. for example to get neutron stars of appropriate mass (see Sec. 
IVD3 of [32]), at each step, we change the value of at the center, in order to make a closer to one. The value of 
4>i (r = 0) to which one wants to converge is modified according to : 

(r = 0) ^ (r = 0) 

Doing so, the value of at the center converges to the only value such that the scaling parameter a is one, providing 
us with the real, non-trivial, solution of the system (17). 

E. Influence of the phcises 

In order to verify that the simplex solver converges to the proper minimum of [^42], we show, on Fig. 8 how various 
quantities vary when one changes the precision at which the extremum is found. Fig. 8 presents the values of : 

• the phase of the second mode ^2 at which the minimum is found, 

• the values of the first modes at the origin 0„ (r = 0), 

• and the value of the minimum of |A2|. 

Those three quantities are shown as a function of w, for three different values of the precision required for the simplex 
solver, 10"-\ 10"'* and 10"^ 

From the values of the phase (p2, it appears that a precision of 10""^ or 10~"* is not good enough, the curve being 
quite noisy. The curve obtained with a precision of 10~^ is smoother and we will assume that this level of precision 
is sufficient for the purpose of this paper and this value will be chosen for all the rest of this work. The situation is 
even better for both the values of the modes at the origin and the actual value of the minimum. Indeed, as can be 
seen on Fig. 8, those quantities show almost not dependence on the level of precision. This is simply related to the 
fact that the values of the fields depend very weakly on the values of the phases of the homogeneous solutions in the 
external region. This is true for the phase of the second mode (p2 but even more for the higher order phases <Pn>2- 
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FIG. 8: Convergence of the various results with respect to the precision required when finding the minimum of \A2\. 

The very moderate dependence of A2 on the phases is ilhistrated by Fig. 9 where the value of the amphtude is 
shown, as a function of the phase ip2, for to = 1.37. AU the other phases (pn>2 are fixed to the same value, each of 
them corresponding to one curve on Fig. 9. It is clear that the influence of the phases of the modes n > 2 is very 
weak. The extrema on the three curves of Fig. 9 are very close to the real extremum found by the simplex solver. 
Therefore we are quite confident that we can find the minimum of \A2\ with a good accuracy, even if the associated 
values of the phases are slightly less accurate. 

F. Convergence tests 

In order to further check the validity and accuracy of our code with respect to the computational parameters, we 

present in this section various additional convergence tests. The idea is to take a reference set of computational 
parameters and to change one of those at a time, in order to check that the obtained results do not change much. 
The first three radial domains are chosen with boundaries [0, 1], [1, 2] and [2,4]. After r = 4, we keep the size of the 
domain constant to 4 (i.e. the boundaries of the domain i arc [4 (« — 3) , 4 (i — 2)]). We foimd that this is necessary 
to ensure that in every domain, we have enough collocation points to resolve the oscillatory homogeneous solutions. 

The various computational parameters are the following : Uz is the number of radial domains which relates directly 
to the value of Rum. given the setting of the domain mentioned above. Nr is the number of coefficients in each domain. 
Finally we will also checked the convergence of the results with respect to the finite number of modes we consider in 
the system (17). 

Our standard setting consists of = 14 => Rum. = 44, Nr = 33, and the use of 6 modes. We then change 
one parameter at a time to verify that the obtained results are indeed meaningful. The convergence is presented by 
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FIG. 9: Value of as a function of the phase ^2, for uj = 1.37. The other phases are maintained equal and set to 1.0 (solid 
line), 1.6 (dashed line) and 2.2 (dotted line). 

showing the same three quantities as in Sec. IV E : the phase ip2, the values of the first modes at the origin and the 

amplitude \A2\. 

The three plots of Fig. 10 show the dependence of the results when varying the number of coefficients in every 
domain. Since for the values of Nr the curves of Fig. 10 change only slightly, this is a strong sign that Nr = 33 is 
greatly sufficient to get accurate results. It is also supported by the fact that, after the end of the iterative scheme 
we find that the relative difference between the right and left hand sides of Eqs. (17) is smaller than 10~®. 

A strong test is provided by the plots of Fig. 11. They show that the results are independent of the value of the 
outer computational radius Rum. thus validating the matching procedure with the oscillatory homogeneous solutions. 
Finally Fig. 12 illustrates the fact that the results do not change substantially when increasing the number of modes. 

The last test of our code has been a comparison of the r.h.s. and l.h.s. of Eq. (3), using the Fourier expansion (16). 
We have computed the relative difference between the r.h.s. and l.h.s. after averaging over one period. The results for 
various number of modes are depicted on Fig. 13 as a function of w. As expected, the error decreases as the number 
of modes increases. This is not surprising, given that Eq. (17) could be satisfied only for an infinite number of modes. 
The error also increases when uj decreases. This can be understood by recalling that the modes are more important 
when CO is small (i.e. see the behavior of (j)n (r = 0) as a function of iv on Figs. 10, 11 and 12). Thus, the effect of the 
missing modes is more important for smaller co. 

Finally we would like to emphasize that the computational parameters exhibited in this section are sufficient to 
compute solutions in the regime of "moderate" frequencies (i.e. frequencies around 1.37 — 1.38). If one wishes to go 
to much lower frequencies, one would need to include more modes, for higher order modes will be more important, 
as indicated by Fig. 13. On the other hand, as will be seen in the next section, the matching point Rum. must be 
increased when one approaches the critical value Uc = \/2. 

V. RESULTS 

As illustrated on Figures 10, 11 and 12 we do find time-periodic quasi-breather solutions of Eqs. (17) for any value 
of the pulsation frequency 1.32 < uj < 1.41, and there is little doubt that such solutions exist for all frequencies in the 
range ]0, It is also clear from the very same Figures that for uj uJc — the family of solutions we consider 
converges pointwise to the trivial solution. We do not expect such quasi-breather solutions of the system (17) to exist 
for u> > LOc- The critical value Uc = is expected to be related to the change of nature of the Helmholtz operator for 

, a/2 being the value at which (pi ceases to decay like an exponential. In the 1 + 1 dimensional case Coron [33] has 
proved that the allowed frequencies are indeed constrained by w < tOc- We have not attempted to prove the analogous 
statement for our case, as it would certainly require some sophisticated mathematical tools. 

The quasi-breather solutions obtained this way are not well localized in space because of their slowly decaying 
cx l/r oscillatory tail. Consequently none of these solutions has finite energy. In fact we have selected a special class 
of time periodic solutions by minimizing the amplitude of their oscillatory tail. This is, in some sense, the closest one 
can get to a breather and we have called those configurations "quasi-breathers" . There is no special value of u) for 
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FIG. 10: Same as Fig. 8 but when varying the Nr in each domain. 



which the aniphtudc of the oscillatory tail would show any tendency to become very small or going to zero. This is 
illustrated on Fig. 14 where the logarithm of the coefhcients of "tails" are shown for the modes 02, (f'Si 4'i s-nd (/)5, as 
a function of uj. The only value of w for which the curves go to zero is the critical one for which the solution tends 
to the trivial one. The many convergence tests of Sec. IV F show that we can accurately compute the value of those 
coefficient and that they are not numerical artifacts. 

To illustrate the behavior of the quasi- breathers we show on Fig. 15 and 16 the energy density ^ and the one 
including the volume element, i.e. r^^ as a function of the radius, for values of w going from 1.36 to 1.40. We can 
clearly see two qualitatively different behaviors: i) the solutions have a well defined "core", where the behavior is 
dominated by the exponential decay of the fields and so the density goes to zero. This core is getting larger when 
U! \f2. ii) However, inevitably at some point, the oscillatory tails start to dominate and the density reaches a 
plateau, ultimately causing the total (integrated) energy to be infinite. Let us mention that if for high values of w, 
the plateau is not seen, it comes solely from the fact that the value of Rxim used on Fig. 15 and 16 is not sufficient. In 
spite of the numerical smallness of the amplitude of this plateau, we are quite confident that its value is reasonably 
accurate, and also that it really corresponds to a physical eff'ect. Indeed, its value is very stable when varying the 
various computational parameters. 

As w increases, the value of the plateau diminishes. This is related to the fact that the coefficients of the homogeneous 
solutions are decreasing functions of a; and to the quadratic nature of the energy density in i.e. the value of the 
plateau is roughly the square of the greatest coefficient (i.e. the one for ^-2). 

As already stated, we can also confirm that the transition radius between the core and the plateau increases as one 
approaches w = \/2. To be more quantitative, we define the transition radius i?trans as the first value of r for which 
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FIG. 11: Same as Fig. 8 but changing the value of Rum.- 
(^1 reaches the amplitude of the oscillations : 

(l)lir = R,n.ns) = —. (22) 

r 

It is the radius at which the oscillatory behavior starts to dominate. iJtrans as a function of w is shown on the left 
panel of Fig. 17. The precise location of the transition radius can not be computed for ui > 1.39, given our nominal 
choice for i?iiin = 44. Once again, this illustrates the fact that when going to higher values of w, one needs to increase 
the boundary radius Rum- On the right panel of Fig. 17 we show the total energy inside the transition radius, i.e. the 
total energy of the core of the configuration. Contrary to the transition radius, this is not a monotonic function and 
there is a configuration of minimum energy. This non trivial behavior is due to two competing effects when increasing 
w, first the increase of the transition radius, and second the decrease of the magnitude of the modes. We also note 
that this minimum is very close to the value of the frequency which Honda and Choptuik [16] claim to be the one 
of their conjectured breather solution (i.e. they quoted ui w 1.366) but this may very well be just a coincidence. It 
seems likely that the energy contained in the core diverges when uj \f2. 

Finally, on Fig. 18 we show the first modes near the origin (left panels) and in the region of the transition radius 
(right panels). Three difi'erent values of w are shown, corresponding to a low value (1.32), the value corresponding to 
the minimum of energy (1.365, cf Fig. 17) and one rather high frequency (1.39). 
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FIG. 13: Relative difference between the l.h.s. and r.h.s. members of Eq. (3). The result is averaged over one period and we 
present the error as a function of w, for 5, 6 and 7 modes. 
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FIG. 15: Energy density for various values of uj, near the origin (left panel) and in the region where the oscillatory tails dominate 
(right panel). i?iim. is equal to 44. 
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FIG. 17: The transition radius between the exponential decay and the oscillatory tail (left panel) and the energy inside this 
radius (right panel), as a function of u). 
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FIG. 18: The first modes near the origin and in the transition region for oj = 1.32, 1.365 and 1.39. 
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VI. THE QUASI-BREATHER CONTENT OF THE OSCILLONS 
A. Fourier decomposition of evolution results 

In this Section we present our oscillon scenario, based on the existence of the quasi-breathers. In particular we 
present some convincing evidence, that an oscillon state of main frequency w is quantitatively described by a quasi- 

brcathcr of the same frequency, whoso tail is cut off at some value of r. In order to do this wc perform Fourier 
decomposition of various a long-lived oscillon solutions found by doing the evolution (see Subsection HE for the 
actual method). This gives us the fundamental pulsation frequency u) and the various modes. Then, those modes 
arc compared to the ones found by directly solving the system (17) for this particular value of lu. The same kind 
of comparison is done in Ref. [16] as to support the existence of a periodic solution of frequency u> « 1.366. An 
important difference is that in Ref. [16] the comparison is performed only for higher frequency near-periodic states in 
an interval of [0 < r < 10] which is comparable to the core part of the quasi-breather and therefore its oscillatory tail 
is not yet apparent. 

It is apparent from the results provided by our evolution code that the solutions corresponding to the time evolution 
data are very close to being time symmetric not only at the center but also in the intermediately far region, justifying 
the choice of only cosine terms in the mode decomposition (16) when constructing QBs. 

B. Near-periodic states 

First we perform the Fourier decomposition of various near-periodic states obtained by fine-tuning the initial 
parameter ro, since these states appear to be periodic and time symmetric to a very high degree. For a specific case 
Figures 19 and 20 show the real part of the modes coming from the Fourier transform of the evolution and the modes 
obtained by the direct solution of (17). The comparison is done with a near-periodic state corresponding to the first 
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FIG. 19: Comparison between the modes coming from solving the system (17) directly and the real part of the modes obtained 
by performing a Fourier transform of the evolution for w = 1.398665. In order to be able to represent on the same graph the 
large values in the central oscillon region along with the small tails at high radii the absolute values are plotted logarithmically. 
We note that the central value (/>n(0) is positive for n = and for odd n while it is negative for other n. The peaks pointing 
downwards on the figure correspond to places where the functions change signature. 

peak on the lifetime curve with initial data (j)c = The Fourier decomposition is performed between the two maxima 
at ti = 1592.29 and t2 = 1596.78. The frequency calculated from the position of the maxima is = 1.398665, which 
is used to calculate the corresponding periodic solution by solving (17). The agreement for the first oscillating mode 
n = 2 is remarkably good even in the oscillating tail region. The curves obtained by the two different decompositions 




24 



0} 



Mode decomposition 
Real part of Fourier transform 




FIG. 20: Comparison of the higher modes of the two systems presented on the previous figure. 



are indistinguishable in the central region. In order to show the excellence of the agreement the radial dependence 
of the relative difference of the corresponding modes obtained by the two methods is shown on Fig. 21. We have 




FIG. 21: Dependence on the radius r of the relative difference of the modes coming from the mode decomposition and the 
Fourier transformation for u = 1.398665. The relative difference is defined as the absolute value of the ratio of the difference 
and the mode decomposition value. 



obtained similarly good agreement not only for u) = 1.398665 but for all frequencies for which near-periodic states 
exist (i.e. for 1.365 < a; < -\/2). 

This comparison between the results coming from the evolution code of Sec. II and from the Fourier-mode decompo- 
sition of the quasi-breathers is a very strong consistency check. Let us recall that the two codes have been developed 
completely independently. 
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Figure 22 shows, for large radii, the first mode for which the oscillatory term appears, i.e. (j)2, in the case lu = 1.38. 
It can be seen from the first of the three plots that the real part of the Fourier transform agrees very well with the 
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FIG. 22: The first two graphs compare the oscillatory region of the mode (f>2 coming from solving the system (17) directly and 
the complex 02 obtained by performing Fourier transform on the result of the evolution code, for uj = 1.38. The first graph 
shows the real and imaginary parts separately, while the second compare the absolute values. The third graph plots the radial 
dependence of the complex argument of the modes 4>i and <p2 obtained by the Fourier transform. 



decomposition of the corresponding QB, the two curves can hardly be distinguished. However, just at those radii 
where the oscillatory behavior appears the imaginary part starts being comparable to the real part, breaking the time 
symmetry in this outer region. The second plot shows that the absolute value of the complex (j)2 obtained by the 
Fourier decomposition essentially behaves like a smooth envelope curve covering the standing wave like peaks of |(^2| 
obtained by the mode decomposition. The presence of the complex argument of the Fourier transform shown on the 
third plot indicates that the oscillating tail obtained by the evolution code is composed of outgoing waves carrying 
out energy from the core region. 



C. Oscillons from Gaussian initial data 



Oscillons obtained by the evolution of generic Gaussian initial data are relatively long living for a large set of the 
possible values of the initial parameter rg. On Fig. 23 wc exhibit the time evolution of the frequency of the two 
typical non fine-tuned states presented on Fig. 1. Although these states are clearly non-periodic, their lifetime is 
still very large compared to the period of the basic high frequency oscillation mode. We have performed the Fourier 
decomposition of these states at various moments of time and compared the results to the mode decomposition of 
the corresponding QB. On Figs. 24, 25, 26 and 27 the Fourier decomposition of the two oscillon states shown on 
Figs. 1 and 23 is given for moments of time where the frequency is approximately 1.30 and 1.36. For comparison the 
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FIG. 23: Time evolution of the underlying high frequency oscillation of the oscillon state of Fig. 1 is presented, 
corresponding modes calculated by solving (17) using the given frequencies is also shown on the figures. Although 
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FIG. 24: The modes of the periodic solution with frequency 1.3 is compared with the Fourier decomposition of two different 
oscillon states. The first oscillon state with initial data 4>c = 1, = 2.70716 is decomposed between the two subsequent 
maxima after t = 411 AT, where the calculated frequency is 1.2995. The second state with cj)^ = \, ro = 2.39 is decomposed just 
after t = 298.72, where the frequency is 1.3011. 

the agreement of the modes is not as excellent as for the near-periodic states, the curves calculated by the two methods 
are still remarkably close. Especially important is the similarity at the tail section of the first oscillating mode, i.e. 
mode 2. This indicates that the periodic quasi-breather solutions can describe, at least qualitatively, even the generic 
long living oscillons. In particular, the existence of the oscillating tail region is responsible for the slow but steady 
energy loss of these configurations. 
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FIG. 25: Comparison of the higher modes of the two systems with frequency uj = 1.3 presented on the previous figure. 
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FIG. 26: The modes of the periodic solution with frequency 1.36 is compared with the Fourier decomposition of two different 
oscillon states. The first oscillon state with initial data 4>c = fo = 2.70716 is decomposed between the two subsequent 
maxima after t = 4158.76, where the calculated frequency is 1.3596. The second state with 0c = 1, = 2.39 is decomposed 
just after t = 1526.58, where the frequency is 1.3613. 



D. Initial data obtained by mode decomposition 

We have also tested our scenario in the reverse way. First, for a given lu (1.38 in the first example), we have solved 
the system (17) for the modes of the quasi-breather. In particular, we can compute (j){t = 0) at the moment of time 
symmetry, and use this as initial data for the evolution code. The evolution of such initial data is shown on Fig. 
28 for numerical simulation with various spatial resolutions. The field (here its value at the origin) oscillates at the 
appropriate frequency for a relatively long period (approx. t=300, with about 66 oscillations) before it collapses in a 
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FIG. 27: Comparison of the higher modes of the two systems with frequency oj = 1.36 presented on the previous figure. 

sub or supercritical way very similar to the collapse of the (f)c = —0.4 states shown on Fig. 4. The initial state is so 
close to the ideal configuration that it even depends on the chosen numerical resolution whether in the final stage the 
configuration collapses in a subcritical or supercritical way. Strangely, the longest living state is not the one with the 
highest resolution, which is probably connected to the fact that the initial data still contain small numerical errors. 
Although the lifetimes attainable by this method are very long compared to the dynamical time scale of the linearized 
problem, they are still shorter than the r « 2000 achievable by fine-tuning the first peak of the cf)c = —0.4 initial data. 
Actually, it should not be surprising that the numerically determined initial data corresponding to a QB cannot be 
as long-living as a configuration which is fine-tuned to 32 decimal digits. 

Since the modes are matched to oscillating tails at large radii, the initial data provided by the mode decomposition 
of a QB is valid up to arbitrarily large r. However, because the slow falloff of the oscillating tails, the energy contained 
in balls of radius r diverges as r goes to infinity. Using such configuration as initial data for our evolution code docs 
not cause serious problems, since the physical distance between grid points in our numerical representation increases 
with the distance from the center, and consequently the high frequency tails cannot be represented numerically 
above a certain radius. This provides a CTitoff in the initial configuration and an effective outer boundary during the 
evolution. An advantage of our conformal compactification method is that this outer boundary moves to higher and 
higher radii when increasing the numerical resolution. On Fig. 29 we show the initial part of the upper envelope curve 
for numerical runs with various spatial resolutions. We emphasize that this strong resolution dependence is entirely 
due to the inadequate representation of the infinite energy initial data. When we used initial data with compact 
support or sufficiently fast falloff, such as the Gaussian initial data in (12), then the code remained convergent up 
to much larger time periods (at the order oit = 10000), and the curves with various resolutions agreed to very high 
precision at the initial stage (at around t = 100). 

An important point is that, with this method, we can generate almost periodic oscillons at frequencies which are 
difficult or impossible to attain by fine-tuning initial data of the form (12). On Fig. 30 we show the upper envelope 
of the oscillations at r = in the initial stage of the evolutions when using initial data provided by the mode 
decomposition method with frequencies lo = 1.30 and lo = 1.36. Again, as the numerical resolution is increased, 
the size of the oscillating tail that can be taken into account in the numerical representation of the initial data gets 
larger, and the evolution will remain nearly periodic for longer times. The field can oscillate truly periodically only 
if the energy lost by radiation through the dynamically oscillating tail is balanced out by an incoming radiation 
already present at the tail section of the initial data. If the tail is cut off above some radius r then the non-periodic 
region moves inwards from that radius essentially at the speed of light. Getting into the non-periodic domain, for 
the frequencies w < 1.36 the tail amplitude appears to be large enough to cause a small but significant energy loss, 
a decrease in amplitude, and consequently a steady frequency increase. On Figs. 31 and 32 we show the long-time 
behavior of the amplitude and frequency of the oscillations generated by initial data with frequencies a; = 1.30 and 
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FIG. 28: On the first graph the upper and lower envelopes of the oscillations of the field at the origin are shown as functions of 
time for numerical simulations with various number of spatial grid points. The initial data is generated by solving the system 
(17) for the modes at the pulsation frequency u! = 1.38. On the second graph the time dependence of the oscillation frequency 
is plotted for each simulation. 



Lo = 1.36. Apart from a relatively short stable initial stage, the main characteristic of the evolution is a slow but 
steady decrease of the amplitude accompanied by a simultaneous increase of the frequency up to a point (approx. 
Lu = 1.365) where the configuration quickly decays. This evolution is very similar to the behavior of generic (i.e. not 
fine-tuned) oscillons started from Gaussian initial data of type (12). A typical example of such evolution, with initial 
data corresponding to an tq close to the top of the lifetime curve (but between two peaks) with (/)c = 1, is also shown 
on Figs. 31 and 32 in order to facilitate comparison. It is also rather remarkable that the onset of the rapid decay 
seems to coincide with the configuration which minimizes the energy inside the core (see Fig. 17). 

The longtime evolution of a low frequency initial data provided by the mode decomposition method is very similar 
to the evolution of generic oscillons evolving from Gaussian initial data. An important difference though, is the smaller 
low frequency modulation of the envelope curve and of the time dependence of the frequency. From this one can expect 
that the Fourier decomposition of these states is even closer to that of the quasi-breather state with the corresponding 
frequency. On Figures 33 and 34 we give the Fourier decomposition of the evolution of the 1.3 initial data between 
the two maxima after t = 5257.45, where the calculated frequency is 1.35997. The agreement is markedly better 
than on Figs. 26 and 27, which is most likely due to the much lower amplitude of the low frequency modulations on 
the envelopes of field value at the center. These low frequency modes arc excited to a higher amplitude by Gaussian 
type of initial data, while they appear to be present to a much lower extent in the lower frequency periodic initial 
data (with uj = 1.3 in this case). We remind the reader that the existence of the peaks on the lifetime curve (see Fig. 
4 of [16]) is closely related to these low frequency oscillations. Namely, the peaks separate domains of the initial data 
in ro with different number of low frequency oscillations on the envelope curve. 
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FIG. 29: Initial behavior of the maximum curve of the evolution of the u) = 1.38 initial data with numerical resolutions 
corresponding to various number of spatial grid points. The points on the curves represent the actual maximum points of the 
oscillating field 4> at the origin r = 0. 
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FIG. 30: Initial behavior of the upper envelope of the oscillations at r = of the evolution of the uj = 1.30 and ui = 1.36 initial 
data for numerical simulations with various number of spatial grid points. 
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FIG. 31: Time evolution of the upper envelope curve of the oscillation of at r = generated by three difFcrent initial data. 
The first two initial data is provided with the mode decomposition method, solving the system (17) by choosing the frequencies 
oj = 1.30 and oj = 1.36. The third evolution corresponds to Gaussian-type initial data (12) with 0c = 1 and ro = 2.70716. 
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FIG. 32: Time evolution of the frequency of the oscillation of at r = for the three kind of evolution simulations shown on 
the previous figure. 



VII. CONCLUSION 



In this paper, we have adapted and applied spectral methods implemented in the LORENE library to find time- 
periodic solutions of the spherically symmetric wave equation of (f/^ theory. Our code passed numerous tests and we 
are quite confident that it is sufficiently precise. With our code we find that for frequencies < a; < \/2 there is 
a whole family of standing wave-type solutions with a regular origin having a well localized core and an oscillatory 
tail. Because of the slow decrease of the oscillatory tail the total energy of these solutions in a ball of sufficiently 
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large radius, R, is proportional to R cf. Figs. (15,16). Wc have constructed a special class of solutions, called quasi- 
breathers, defined by minimizing the amplitude of the oscillatory tail of the solutions. The size of the core of the QBs 
gets larger and larger as a; — > \/2 and the amplitude of the oscillatory tail is getting smaller and smaller. Interestingly 
the energy contained in the core of the QBs exhibits a minimum for lu « 1.365. 

Using our high precision time evolution code, we have also investigated the time evolution of Gaussian initial data. 
We observe that a generic oscillon state can be characterized by a slowly varying frequency, oj{t), increasing in time; 
up to a critical one Uc, when it decays rapidly. We have also investigated in detail the near-periodic states described 
first in Ref. [16] which can be characterized by an almost constant frequency. In contradistinction to Ref. [16] we find 
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such near-periodic states for any frequency \/2 > cj > ojc- Moreover we find that the near-periodic states decay in 
time (i.e. they cannot be truly time periodic states) cf. Figs. (5,6). In particular while loosing some of their energy 
their frequency decreases very slowly with time. 

By a careful comparison of the Fourier modes of an oscillon state of frequency Lo{t) with that of the corresponding 
QB, we have obtained convincing evidence that the localized part of the oscillon is nothing but the core of the QB of 
the same frequency. Our results demonstrate that the time evolution of an oscillon state can be described to a good 
approximation as an adiabatic evolution through a sequence of QBs with a slowly changing frequency Lo{t). What is 
more, the oscillatory tail of the QB describes very well the standing wave part of the oscillon. Therefore any oscillon 
contains the core and a significant part of the oscillatory tail of the corresponding QB. The existence of near-periodic 
states is closely related to the fact that the energy of the oscillon core exhibits a minimum for a critical frequency 
cOc ~ 1.365 and for \/2 > w > Wc a single unstable mode appears which can be suppressed by fine-tuning the initial 
data. 
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